Cite this Article
Hou Zhang-lin, Ju Ying, Zong Yi-wu, Ye Fang-fu, Zhao Kun. Molecular dynamics simulations on the dynamics of two-dimensional rounded squares. Chinese Physics B, 2018, 27(8): 088203  
Permissions
Molecular dynamics simulations on the dynamics of two-dimensional rounded squares
Hou Zhang-lin1, Ju Ying1, Zong Yi-wu1, Ye Fang-fu2, ‡, Zhao Kun1, §
Key Laboratory of Systems Bioengineering (Ministry of Education), School of Chemical Engineering and Technology, Tianjin University, Tianjin 300072, China
Beijing National Laboratory for Condensed Matter Physics and Key Laboratory of Soft Matter Physics of Chinese Academy of Sciences, Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China

 

† Corresponding author. E-mail: fye@iphy.ac.cn kunzhao@tju.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 21573159 and 21621004).

Abstract

The collective motion of rounded squares with different corner-roundness ζ is studied by molecular dynamics (MD) simulation in this work. Three types of translational collective motion pattern are observed, including gliding, hopping and a mixture of gliding and hopping. Quantitatively, the dynamics of each observed ordered phase is characterized by both mean square displacement and van Hove functions for both translation and rotation. The effect of corner-roundness on the dynamics is further studied by comparing the dynamics of the rhombic crystal phases formed by different corner-rounded particles at a same surface fraction. The results show that as ζ increases from 0.286 to 0.667, the translational collective motion of particles changes from a gliding-dominant pattern to a hopping-dominant pattern, whereas the rotational motion pattern is hopping-like and does not change in its type, but the rotational hopping becomes much more frequent as ζ increases (i.e., as particles become more rounded). A simple geometrical model is proposed to explain the trend of gliding motion observed in MD simulations.

PACS: 82.70.Dd;64.75.Yz
Keyword:molecular dynamics simulation;rounded square;hopping;gliding;collective motion
1. Introduction

Phase behavior of two-dimensional (2D) colloidal systems including both static structures and dynamics depend in an important way on the geometrical shape of colloids. For hard particles with a shape as simple as a disk, besides solid and liquid phases, they can also form a liquid crystalline hexatic phase, which has quasi-long-range bond orientational order but short-range positional order.[1–4] The Kosterlitz–Thouless–Halperin–Nelson–Young (KTHNY) theory predicts a two-step melting scenario from solid to liquid via the hexatic phase.[2–4] Recent simulation work by Bernard and Krauth[5] confirmed the two-step melting process and found that the solid-to-hexatic transition is continuous while hexatic-to-liquid transition is first-order.

Dynamic properties of 2D hard colloidal systems have also been studied extensively. In a dense quasi-2D liquid state, Marcus et al.[6] showed a hopping motion of particles, which can cause the single-particle displacement distribution to deviate from Gaussian form. Zangi and Rice[7] further studied the cooperative motion in such systems with molecular dynamics (MD) method, and found the correlation between the hopping motion and dynamic heterogeneity. They showed that the deviation from Gaussian form starts from liquid phase and persists all the way through hexatic phase into solid phase. Kim et al.[8] found that 2D hard discs in solid and hexatic phases have seemingly Fickian (i.e., the mean-square displacement is linear with time) but strongly heterogeneous dynamics indicated by the non-Gaussian behavior of the self-part of the van Hove correlation function. They showed the dynamic heterogeneity is closely related to the cage formation/breakdown and the hopping motion of particles. Zahn and Maret[9] studied the dynamics of a quasi-2D system of superparamagnetic colloidal particles at the melting transition, and found that the dynamics in hexatic phase is governed by the glide of free dislocations.

Compared with spheres or disk-like particles, colloids with non-spherical (or non-discotic) shapes, such as triangles, squares, pentagons, and crosses show much richer phase behavior in two dimensions.[1,10–23] For example, a new liquid crystalline triatic phase has been found in a 2D triangle system; two square crystals with opposite chirality were also observed in a 2D square cross system. To systematically study the effect of particle shape on the phase behavior of anisotropic particles, Anderson et al.[1] performed a Monte Carlo (MC) simulation for a series of hard regular polygons including triangles all the way up to 14-gons, and displayed a rich phase diagram of particles. They showed that the melting transition depends upon both particle shape and symmetry considerations.

The phase behavior of anisotropic particles is very sensitive to the details of particle shape. For a 2D square system, the experimental results showed a phase sequence of liquid-rotator hexagonal crystal-rhombic crystal,[11] which is completely different from the liquid-tetratic-square crystal obtained in a MC simulation of a hard square system.[17,20] The discrepancy is due to the corner rounding of squares used in the experiment, which was confirmed in a detailed MC simulation of rounded hard squares by Avendaño and Escobedo.[20]

Most of previous studies on the non-spherical (or non-discotic) colloidal systems have been focused on the static structures. The dynamics of 2D anisotropic colloidal systems, however, are much less studied. Particularly, how particle shape affects the dynamics of 2D colloidal systems consisting of anisotropic particles remains largely unexplored.

In this work, we report a systematic study on dynamics of rounded squares with different corner-roundness that interact through a short-range repulsive interaction using MD simulations. For systems in which particles interact with hard particle interactions or short-range repulsions, considering that the effect of particle shape will be more prominent when particles are in close proximity to each other, we focus on the dynamics of ordered phases (including liquid crystalline phases) which are at relatively large surface fractions (defined as the ratio of the surface area occupied by all particles in a simulation box to the total area of the simulation box). Mean square displacements and van Hove functions are calculated to characterize the dynamics of each observed ordered phase. The effect of corner-roundness on the dynamics is further studied by comparing the dynamics of rhombic crystal phases formed by different corner-rounded particles at a same surface fraction. A simple model is also proposed to explain the trend of gliding motion observed in MD simulations.

This paper is organized as the following: in Section 2, the model and simulation methods used in this work are introduced; in Section 3, results are shown and discussed; and finally, a conclusion is given in Section 4.

2. Model and simulation methods

A coarse-grained model of rounded squares is used in this work. First, the outline of a rounded square is consisting of four straight line segments connected with four arcs. Each arc is obtained by cutting a corner of a perfect square (of length L) with a circle of radius R0 tangent to the two associated sides of the corner, as illustrated in the left panel of Fig. 1(a) (see Supplementary Information for more details). The degree of roundness of a particle is defined as ζ = R0/Rmax, where Rmax = L/2 is the radius of the inscribed circle of the square. So ζ = 0 corresponds to a square shape, and ζ = 1 corresponds to a disk shape. To simplify the particle interaction calculation, we use a coarse-grained model to represent the rounded square by putting spheres with a diameter σ (σ = L/23) along the sides of the rounded square. The spheres are fixed with the rounded square and equally distributed with an interval separation distance of 0.5σ along the sides, and interact with a Weeks–Chandler–Andersen (WCA) potential,[24,25] which is the 12-6 Lennard-Jones potential truncated at the minimum rij = 21/6σ and shifted vertically by ε to be purely repulsive,[25] where rij is the distance between two spheres, ε is the characteristic energy. Using this coarse-grained model, the interaction between two rounded squares is calculated by summing over the pair potential of each possible pair of spheres between the two particles (pairs of spheres from the same rounded square are excluded).

Fig. 1. (color online) (a) Sketch of construction of a corner-rounded square (left), and a coarse-grained model of the rounded square used in MD simulations (right). (b) Phase diagram of rounded squares in the η–ζ plane obtained by MD simulations with the NVT ensemble. I: isotropic liquid phase (open circle); H: hexatic phase (open star); RHX: hexagonal rotator crystal phase (solid triangle); T: tetratic phase (open square); RB: rhombic crystal (open diamond).

MD simulations of rounded squares are performed in the NVT ensemble and Nose-Hoover thermostat is employed to maintain the temperature T of simulated systems[26–31] using a GPU-based MD package, the GPU-Accelerated Large-scale Molecular Simulation Toolkit (GALAMOST).[32,33] The reported length, energy, temperature T, and time t are in units of σ, ε, ε/kB, and , respectively. m is the particle mass. A lower T corresponds to a larger ε relative to kBT, and at a lower T the pair interaction of particles would be more closer to the hard particle interaction used in MC simulations[20] and experiments.[11] However, it also takes a longer time for the system to get equilibrated at a lower T. So, based on those considerations, MD simulations are performed at T = 0.1 in this work. For each simulation, periodic boundary conditions are applied and the initial configuration is taken to be an appropriate rhombic lattice (see Table S1 in Supplementary Information) with a total number of 3264 ∼ 4600 particles. The time step is Δt = 0.002. First, the system run at least 107 MD steps at a higher temperature T = 1.0 to let particles on the initial rhombic lattice quickly diffuse and be uniformly distributed in the simulation box. The resulted configuration of particles at the end of each run at T = 1.0 would serve as the starting point for subsequent runs at T = 0.1. Then, the system run another 0.5 × 107 MD steps at T = 0.1 to equilibrate. After that, the system continued to run 2.0 × 107 MD steps at T = 0.1 to output configurations every 5 × 104 MD steps for statistical analysis of dynamics (for static structures, the system just run 0.1 × 107 MD steps). In this work, the simulation results for rounded squares with ζ = 0.09, 0.286, 0.4, 0.5, 0.667, and 0.8 are shown.

To characterize the diffusive behavior of particles, the mean square displacement (MSD) MSD ≡ ⟨(Δ ri(t))2⟩ and the mean square angular displacement (MSAD) MSAD ≡ ⟨(Δ θi(t))2⟩ are calculated from the tracked trajectories of particles. Here Δri(t) = ri(t) − ri (0) and Δθi (t) = θi (t) − θi (0) are relative displacements of particle i during time interval t for translation and rotation, respectively. ⟨···⟩ denotes ensemble average.

Another useful parameter to characterize the dynamics of rounded squares is the self-part of the van Hove correlation function, , which gives the probability that the displacement Δri(t) is equal to distance r at time interval t. In addition, the distinct-part of the van Hove correlation function, with Δrij(t) = rj(t) − ri (0)[34] is also calculated, which describes the average motion of the remaining N − 1 particles (excluding the reference particle i) relative to the position of particle i at t = 0. Gd(r, t) is related to the pair correlation function g(r) by Gd(r, 0) = ρg(r), where ρ is the number density. Similarly for rotational motion, the angular van Hove correlation function is defined as .

To better characterize the dynamics of 2D solids, an alternative way to measure the self-part of the van Hove correlation function is to replace Δri(t) by the relative displacement of neighboring particles (indices i and i + 1) Δrrel(t) = Δ ri + 1 (t) − Δri (t),[9,35] to get

3. Results and discussion
3.1. The phase diagram obtained using MD and collective motion of particles

The phase diagram obtained by MD simulations is shown in Fig. 1(b). There are five phases observed, which are isotropic (I), hexatic (H), tetratic (T), hexagonal rotator crystal (RHX), and rhombic crystal (RB) phases. The detailed results of the order parameters and the correlation functions associated with each phase are in the Supplementary Information (Figs. S1–S7). These phases are also observed using MC simulations,[20] although more MC simulations with larger systems are needed for a conclusive result of hexatic phase. However, the MD results and the MC results also show some differences for rounded squares with specific roundness, particularly for ζ = 0.09, 0.286, and 0.4. For the systems of ζ = 0.09 and ζ = 0.286, a square crystal phase is observed in MC simulations but not in MD simulations. For the system of ζ = 0.4, the MD results show a phase sequence of I–H–RB while the MC results display a sequence of I–RHX–RB. Those discrepancies might be originated from the different models of particles employed in MC and MD simulations. In MC simulations, the particles are mathematically ideal rounded squares with hard particle interactions while in MD simulations, a coarse-grained model of rounded squares is used and particles interact with a WCA-based pair potential. But a detailed analysis on the differences between MC and MD simulations is beyond the scope of this work.

Dynamic analysis in the observed ordered phases at high surface fractions reveals three types of translational collective motion pattern including hopping, gliding and a mixture of hopping and gliding, as shown in Fig. 2. Hopping motion pattern is observed typically in ordered phases of rounded squares with high ζ, such as in RHX and RB phases of ζ = 0.667, and in H and RHX phases of ζ = 0.8. Figure 2(a)–2(d) show an example of hopping motion at different time scales in the RHX phase of ζ = 0.8 at surface fraction η = 0.73. At t = 0.5 × 104, shows a single peak at r/L < 1. This means that, at this time scale, most particles move within the cages formed by their neighbors, and only a few particles show a maximum displacement larger than L, as shown in Fig. 2(b). At t = 2.5 × 104, more particles show a larger displacement (Fig. 2(c)) and a second peak (at 1 < r/L < 2) of begins to develop, which becomes more apparent at t = 3.5 × 104, indicating an increased hopping probability. Consequently, at t = 3.5 × 104 the hopping clusters grow and can merge into large ones (Fig. 2(d)). In addition, Pd (r,t) ≡ Gd(r,t)/ρ shows multiple peaks at roughly the same peak locations for all tested time scales, indicating that the structure of the simulated system does not change much. More interestingly, at the long time scale, the location of the minimum in the region of r/L < 1 is not at r = 0 (indicated by the olive arrow in the inset of Fig. 2(a)). This indicates that the probability of finding a particle at a location where another particle resides at t = 0 increases,[7] and it is higher than the probability to find particles at places around that location. By following particles motion in a hopping process, we found that particles perform a hopping motion with the help of transient defects formed around the hopping particles (Movie S1).

Fig. 2. (color online) Three types of translational collective motion pattern. (a) An example of hopping motion pattern (RHX phase of ζ = 0.8 at η = 0.73) characterized by and Pd(r,t) ≡ Gd(r,t)/ρ. An example of a hopping trajectory is shown in black circle inset (different color represents different particle); (b)–(d) Snapshots of a representative trajectory at 3 different times t = 0.5 × 104 (b), 2.5 × 104 (c), and 3.5 × 104 (d). (e) An example of gliding motion pattern (T phase of ζ = 0.09 at η = 0.79) characterized by Ps(r,t) ≡ 2πrGs (r,t) and Pd(r,t). An example of a gliding motion is shown in black circle inset (different color represents different particle); (f)–(h) Snapshots of a representative trajectory at t = 0.3 × 104 (f), 1.0 × 104 (g), and 3.0 × 104 (h). (i) An example of a mixture of hopping and gliding motion pattern (H phase of ζ = 0.5 at η = 0.73) characterized by Ps(r,t) and Pd(r,t). (j)–(l) Snapshots of a representative trajectory at t = 0.5 × 104 (j), 2.5 × 104 (k), and 3.5 × 104 (l). The color code of each trajectory is determined based on the maximum displacement of the corresponding particle in a given time.

Gliding motion pattern is typically observed in T phase and RB phase with low degree of roundness (ζ ≤ 0.4). Figures 2(e)–2(h) show an example of gliding motion at different time scales in T phase of ζ = 0.09 at η = 0.79. At t = 0.3 × 104, Ps(r,t) ≡ 2πrGs(r,t) shows a single peak. But unlike the hopping motion, trajectories of particles (Fig. 2(f)) show a stripe-like pattern which is the characteristic of gliding. At this time scale, the multiple peaks in Pd(r, t) that exist at t = 0.1 × 104 disappear, and Pd(r,t) at r/L < 1 is elevated but still below 1, indicating that due to gliding of particles, local structures have been changed but the displacement under the current time scale is still not big enough to let all particles completely lose their original positions at t = 0. At t = 1.0 × 104, Ps(r, t) shows a second peak, which is separated from the first peak by 0.8L. The trajectories of particles at this time scale show a grid-like pattern that percolates the simulation box as shown in Fig. 2(g) and Pd(r, t) shows a flat curve implying that the initial configuration is lost completely due to gliding. At t = 3.0 × 104, the first peak of Ps(r, t) shifts to a larger distance indicating the dynamics of the whole system is liquid-like. At this long time scale, the trajectories in Fig. 2(h) show that the gliding clusters grow into one big network which encloses some small islands consisting of particles with relatively small displacements.

In H phases of ζ ≤ 0.667, and crystal phases of ζ = 0.5, a mixture of hopping and gliding motion is observed. Figures 2(i)–2(l) show an example of mixed motion at different time scales in a H phase of ζ = 0.500 at η = 0.730. At t = 0.5 × 104, Ps(r, t) shows a single peak, and the trajectories of particles in Fig. 2(j) show both local compacted clusters and stripe-like clusters indicating a mixture of hopping and gliding. At this time scale, Pd(r, t) shows multiple peaks indicating the initial configuration has not been completely disrupted yet. At t = 1.0 × 104, a second peak develops in Ps(r, t) and becomes even stronger than the first peak at t = 2.5 × 104 and 3.5 × 104, indicating a stronger degree of dynamic heterogeneity. The trajectories at t = 2.5 × 104 and 3.5 × 104 show a broad connected region of particles with maximum displacements larger than L, decorated with a few separated branched regions of particles with maximum displacements larger than 2L. Peaks of Pd(r, t) vanish for t > 1.0 × 104 due to gliding motion, and Pd(r, t) at r/L < 1 increases with time scale and approaches to be flat, indicating the initial configuration is lost as particles glide to large displacements. Note that in the motion patterns involving gliding motion, we use Ps(r, t) rather than . The reason is that in some systems in which particles can have both gliding motion and hopping motion, if is used, the gliding motion will mask the effect of the hopping motion by smearing out relatively small peaks found in Ps(r, t). This can be clearly seen in Fig. S8 that shows the results using from the same data as in Fig. 2(i).

3.2. Effect of corner-roundness on dynamics

To understand the effect of particle shape on dynamics, we compare the motion of particles with different corner-roundness. We choose the RB phases of ζ = 0.286, 0.4, 0.5, and 0.667 which are all at the same surface fraction η = 0.83 in order to minimize the effect of different phases and different surface fractions on dynamics. The lattice angle of RB phases at η = 0.83 is 82.7°, 77.7°, 70.5°, and 68.0° for the systems of ζ = 0.286, 0.4, 0.5, and 0.667, respectively. Results are shown in Fig. 3.

Fig. 3. (color online) Dynamics of RB phases of ζ = 0.286, 0.4, 0.5, and 0.667 at the same surface fraction η = 0.83. (a) MSD of RB phases; (b)–(e) Ps(r,t) and Pd(r,t) for RB phase of ζ = 0.286 (b), 0.4 (c), 0.5 (d), and 0.667 (e). Insets at the bottom right are snapshots of a representative translational trajectory at t = 3.5 × 104. The color code of each trajectory is determined based on the maximum displacement of the corresponding particle in a given time. (f) MSD of RB phases; (g)–(j). Gs(θ,t) for RB phases of ζ = 0.286 (g), 0.4 (h), 0.5 (i), and 0.667 (j). Insets at the upper right are snapshots of RB phases at t = 3.0 × 104, in which each particle is color-coded based on its maximum angular displacement in a given time.

As ζ changes from 0.286 to 0.667, the slope of MSD curve decreases from above 1 to below 1, and for ζ = 0.5 and 0.667, the MSD curves tend to be flat at relatively long time scale (Fig. 3(a)), indicating that the translational motion of particles has changed greatly from a super-diffusive motion to a confined diffusion. By contrast, for rotational motion, as ζ increases the MSAD curves (Fig. 3(f)) are very similar and only the magnitude of MSAD increases, indicating that the rotational dynamics becomes faster as particles become more rounded. Correspondingly, the rotational motion patterns of different rounded squares are also similar and only hopping motion is observed, which contributes to the increase in MSAD. However, their translational motion patterns show a large difference. For systems of ζ = 0.286 and ζ = 0.4, gliding motion is dominant. Trajectories (inset in Figs. 3(b) and 3(c)) exhibit strong stripe-like patterns and particles can migrate a very long distance (most particles can migrate at least 5L in the ζ = 0.286 system, and 3L in the ζ = 0.4 system). The peaks of Pd(r, t) for both systems decrease rapidly with time and disappear at t = 0.4 × 104. The second peak of Ps(r, t) is due to more particles involved in the large gliding clusters, similar to the case in T phase (Fig. 2(b)). But unlike the relatively large displacement of the first peak of Ps(r, t) in T phase (which can shift continuously to almost 1L at t = 3.0 × 104), in RB phase, the first peak of Ps(r, t) fluctuate around the position of 0.3L ∼ 0.4L. This implies that the dynamics of gliding motion in RB phase is slower than that in T phase, which is likely due to the different ordering between RB and T phases. Because of the gliding motion, the MSD curves show a super-diffusive behavior for t < 1.0 × 104. For ζ = 0.5, both hopping and gliding motion are observed (see the snapshot of a representative translational trajectory in the inset of Fig. 3(d)). At t = 2.0 × 104 the peaks in Pd(r, t) disappear and Pd(r, t) at r/L < 1 increases, which are characteristics of the gliding motion. The hopping motion also contributes to the appearance of the second peak in Ps(r, t). For ζ = 0.667 system (Fig. 3(e)), only hopping motion is occasionally observed, which is typically associated with the defects in RB crystals. Most particles are confined and vibrate around their equilibrium positions, which leads to the plateau of the MSD curve. A similar trend of particle motion is also observed when RHX phases or H phases of different roundness are used for comparison (see Figs. S9 and S10 in Supplementary Information).

Those results suggest that the corner-roundness will affect more on translational motion than on rotational motion. As particles become more rounded, translationally the gliding motion is gradually suppressed and the hopping motion becomes dominant, whereas rotationally the rotational hopping becomes more frequent (we can imagine that for a limiting case, hard frictionless disks can freely rotate even at highest packing density). The different behavior of gliding and hopping observed in different rounded particles might have a simple geometrical reason. Considering the ζ = 0.286 particle whose sides still have a relatively large straight portion, when those particles form a RB lattice, the side-side contact (i.e., straight line-straight line contact) between particles will be still dominant, so when particles are aligned in the RB lattice, their sides can form relatively smooth boundaries between rows/columns of the lattice (Fig. 4(a)), which will facilitate the gliding motion. By contrast, for the ζ = 0.667 particles whose sides have a considerable arc potion, when they form a RB lattice, the side-point contact (i.e., arc-straight line contact) becomes more frequent, and more importantly, their sides form very rough boundaries (e.g., with a zigzag shape) between rows/columns of the lattice (Fig. 4(b)), which can thus suppress the gliding motion. This model agrees with the structures obtained in the corresponding simulations (Figs. 4(c) and 4(d)).

Fig. 4. (color online) A geometrical model showing (a) a smooth boundary (green line) for ζ = 0.286 rounded squares and (b) a zigzag boundary (red line) for ζ = 0.667 rounded squares between two rows of a rhombic lattice along one lattice direction at η = 0.830; (c)–(d) Snapshots from simulations of RB phase at η = 0.830 show (c) a smooth boundary (green line) along which gliding is observed for ζ = 0.286 rounded squares and (d) rough boundaries (red lines) and thus no gliding motion for ζ = 0.667 rounded squares.
4. Conclusion

We investigate the dynamics of rounded squares that interact through a WCA potential using MD simulations. Three types of motion pattern, including gliding, hopping and a mixture of gliding and hopping are found. Gliding motion is typically observed in T phase, as well as in RB phases of ζ ≤ 0.4. Whereas, in ordered phases of large ζ, such as in RHX and RB phases of ζ ≥ 0.667, hopping motion pattern is typically observed. The third motion pattern is a mixture of gliding and hopping, which is typically observed in crystal phases of ζ = 0.5, as well as in H phase of ζ ≤ 0.667.

The effect of corner-roundness on dynamics is also studied. The results show that as ζ increases from 0.286 to 0.667 in RB phases, translational motion changes from being super-diffusive with a slope of MSD larger than 1 to a confined diffusion with a MSD curve that can reach to a plateau within the tested time scales. Accordingly the translational collective motion of particles changes from a gliding dominant pattern to a hopping dominant pattern. By contrast, the rotational motion pattern does not change in its type and is hopping-like from ζ = 0.286 to ζ = 0.667, but the rotational hopping becomes more and more frequent as ζ increases. A simple geometrical model can be used to better understand the observed different behavior of gliding and hopping in systems of different ζ. Our work will shed light on the relationship of particle shape, local structures and dynamics.

Reference
[1] Anderson J A Antonaglia J Millan J A Engel M Glotzer S C 2017 Phys. Rev. X 7 021001
[2] Kosterlitz J M Thouless D J 1973 J. Phys. C: Solid State Phys. 6 1181
[3] Halperin B I Nelson D R 1978 Phys. Rev. Lett. 41 121
[4] Young A P 1979 Phys. Rev. B 19 1855
[5] Bernard E P Krauth W 2011 Phys. Rev. Lett. 107 155704
[6] Marcus A H Schofield J Rice S A 1999 Phys. Rev. E 60 5725
[7] Zangi R Rice S A 2004 Phys. Rev. Lett. 92 035502
[8] Kim J Kim C Sung B J 2013 Phys. Rev. Lett. 110 047801
[9] Zahn K Maret G 2000 Phys. Rev. Lett. 85 3656
[10] Zhao K Bruinsma R Mason T G 2012 Nat. Commun. 3 801
[11] Zhao K Bruinsma R Mason T G 2011 Proc. Natl. Acad. Sci. USA 108 2684
[12] Zhao K Mason T G 2009 Phys. Rev. Lett. 103 208302
[13] Zhao K Mason T G 2015 Proc. Natl. Acad. Sci. USA 112 12063
[14] Zhao K Mason T G 2014 J. Phys.: Condens. Matter 26 152101
[15] Zhao K Mason T G 2012 J. Am. Chem. Soc. 134 18125
[16] Rossi L Mason T G 2015 Soft Matter 11 2461
[17] Wojciechowski K W Frenkel D 2004 Comput. Meth. Sci. Technol. 10 235
[18] Schilling T Pronk S Mulder B Frenkel D 2005 Phys. Rev. E 71 036138
[19] Benedict M Maguire J F 2004 Phys. Rev. B 70 3352
[20] Avendaño C Escobedo F A 2012 Soft Matter 8 4675
[21] Zhang T H Cao J S Liang Y et al. 2016 Acta Phys. Sin. 65 176401 in Chinese
[22] Feng Wang Han Y L 2018 Physics 47 238 in Chinese
[23] Sun Y L Wang H G Zhang Z X 2018 Acta Phys. Sin. 67 106401 in Chinese
[24] Weeks J D 1971 J. Chem. Phys. 54 5237
[25] Ahmed A Sadus R J 2009 Phys. Rev. E 80 061101
[26] Nosé S 1984 J. Chem. Phys. 81 511
[27] Nosé S 2002 Mol. Phys. 100 191
[28] Hoover W G 1985 Phys. Rev. A 31 1695
[29] Evans D J Holian B L 1985 J. Chem. Phys. 83 4069
[30] Hoover W G 1986 Phys. Rev. A 34 2499
[31] Frenkel D Smit B 2001 Understanding Molecular Simulation: From Algorithms to Applications New York Academic Press, Inc. 147 158
[32] Zhu Y L Liu H Li Z W Qian H J Milano G Lu Z Y 2013 J. Comput. Chem. 34 2197
[33] Zhu Y L Pan D Li Z W Liu H Qian H J Zhao Y Lu Z Y Sun Z Y 2018 Mol. Phys. 116 1065
[34] Hopkins P Fortini A Archer A J Schmidt M 2010 J. Chem. Phys. 133 249
[35] Han Y Ha N Y Alsayed A M Yodh A G 2008 Phys. Rev. E 77 041406